Wykorzystane biblioteki:
library(ggplot2) #Do tworzenia regularnych wykresów
library(plotly) #Do tworzenia interaktywnych wykresów
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library(dplyr) #Do manipulacji na danych
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(magrittr) #Dla operatora %>%
library(caTools)
library(reshape2) #Dla funkcji 'melt' upraszczającej stworzenie heatmapy korelacji
library(knitr) #Dla prezentacji dataframe'a w sensowny sposób w RMarkdown
Wymuszanie identycznych rezultatów dla kolejnych wywołań i wczytanie danych:
set.seed(2112) #Zadbanie o taki sam rezultat dla kolejnych wywołań
data <- read.csv('Life_Expectancy_Data.csv')
Funkcje pomocniczne do pokazania podsumowań:
sum_nas <- function(x, na.rm){
sum(is.na(x))
}
statistics <- list(list('Nothing_to_see_here', "Column name"),
list(mean, "Mean"), list(sd, "Standard Deviation"),
list(min, "Minimum"), list(max, "Maximum"),
list(median, "Median"), list(sum_nas, "Amount of NaN"))
all_cols = list()
for (x in statistics){
all_cols[[length(all_cols)+1]] <- x[[2]]
}
df_summary <- data.frame(matrix(ncol = length(all_cols), nrow = 0))
colnames(df_summary) <- all_cols
Przedstawienie danych ogólnych o datasecie i zmiennych kategorycznych:
cols <- ncol(data)
cat("Number of rows: ", nrow(data), "\nNumber of columns: ", cols, "\n")
## Number of rows: 2938
## Number of columns: 22
ite <- 1
for (x in data){
dist = " "
if (class(x) != 'factor'){
new_row <- list()
for (y in statistics){
if (y[[2]] == 'Column name')
new_row[[length(new_row)+1]] <- colnames(data)[ite]
else
new_row[[length(new_row)+1]] <- y[[1]](x, na.rm=TRUE)
}
df_summary[nrow(df_summary) + 1, ] <- new_row
}
else{
names <- unique(x)
cat('\nFactor column: ', colnames(data)[ite], ', Unique values: ', paste(shQuote(names, type="cmd"), collapse=", "), '\n')
}
ite <- ite+1
}
##
## Factor column: Country , Unique values: "Afghanistan", "Albania", "Algeria", "Angola", "Antigua and Barbuda", "Argentina", "Armenia", "Australia", "Austria", "Azerbaijan", "Bahamas", "Bahrain", "Bangladesh", "Barbados", "Belarus", "Belgium", "Belize", "Benin", "Bhutan", "Bolivia (Plurinational State of)", "Bosnia and Herzegovina", "Botswana", "Brazil", "Brunei Darussalam", "Bulgaria", "Burkina Faso", "Burundi", "Cote d'Ivoire", "Cabo Verde", "Cambodia", "Cameroon", "Canada", "Central African Republic", "Chad", "Chile", "China", "Colombia", "Comoros", "Congo", "Cook Islands", "Costa Rica", "Croatia", "Cuba", "Cyprus", "Czechia", "Democratic People's Republic of Korea", "Democratic Republic of the Congo", "Denmark", "Djibouti", "Dominica", "Dominican Republic", "Ecuador", "Egypt", "El Salvador", "Equatorial Guinea", "Eritrea", "Estonia", "Ethiopia", "Fiji", "Finland", "France", "Gabon", "Gambia", "Georgia", "Germany", "Ghana", "Greece", "Grenada", "Guatemala", "Guinea", "Guinea-Bissau", "Guyana", "Haiti", "Honduras", "Hungary", "Iceland", "India", "Indonesia", "Iran (Islamic Republic of)", "Iraq", "Ireland", "Israel", "Italy", "Jamaica", "Japan", "Jordan", "Kazakhstan", "Kenya", "Kiribati", "Kuwait", "Kyrgyzstan", "Lao People's Democratic Republic", "Latvia", "Lebanon", "Lesotho", "Liberia", "Libya", "Lithuania", "Luxembourg", "Madagascar", "Malawi", "Malaysia", "Maldives", "Mali", "Malta", "Marshall Islands", "Mauritania", "Mauritius", "Mexico", "Micronesia (Federated States of)", "Monaco", "Mongolia", "Montenegro", "Morocco", "Mozambique", "Myanmar", "Namibia", "Nauru", "Nepal", "Netherlands", "New Zealand", "Nicaragua", "Niger", "Nigeria", "Niue", "Norway", "Oman", "Pakistan", "Palau", "Panama", "Papua New Guinea", "Paraguay", "Peru", "Philippines", "Poland", "Portugal", "Qatar", "Republic of Korea", "Republic of Moldova", "Romania", "Russian Federation", "Rwanda", "Saint Kitts and Nevis", "Saint Lucia", "Saint Vincent and the Grenadines", "Samoa", "San Marino", "Sao Tome and Principe", "Saudi Arabia", "Senegal", "Serbia", "Seychelles", "Sierra Leone", "Singapore", "Slovakia", "Slovenia", "Solomon Islands", "Somalia", "South Africa", "South Sudan", "Spain", "Sri Lanka", "Sudan", "Suriname", "Swaziland", "Sweden", "Switzerland", "Syrian Arab Republic", "Tajikistan", "Thailand", "The former Yugoslav republic of Macedonia", "Timor-Leste", "Togo", "Tonga", "Trinidad and Tobago", "Tunisia", "Turkey", "Turkmenistan", "Tuvalu", "Uganda", "Ukraine", "United Arab Emirates", "United Kingdom of Great Britain and Northern Ireland", "United Republic of Tanzania", "United States of America", "Uruguay", "Uzbekistan", "Vanuatu", "Venezuela (Bolivarian Republic of)", "Viet Nam", "Yemen", "Zambia", "Zimbabwe"
##
## Factor column: Status , Unique values: "Developing", "Developed"
Przedstawienie informacji o zmiennych numerycznych dataseta:
#Wzorzec: https://www.rdocumentation.org/packages/knitr/versions/1.36/topics/knit_print
knit_print.data.frame = function(x, ...) {
res = paste(c("", "", kable(x, output = FALSE)), collapse = "\n")
asis_output(res)
}
knit_print(df_summary)
| Column name | Mean | Standard Deviation | Minimum | Maximum | Median | Amount of NaN |
|---|---|---|---|---|---|---|
| Year | 2.007519e+03 | 4.613841e+00 | 2000.00000 | 2.015000e+03 | 2.008000e+03 | 0 |
| Life.expectancy | 6.922493e+01 | 9.523867e+00 | 36.30000 | 8.900000e+01 | 7.210000e+01 | 10 |
| Adult.Mortality | 1.647964e+02 | 1.242921e+02 | 1.00000 | 7.230000e+02 | 1.440000e+02 | 10 |
| infant.deaths | 3.030395e+01 | 1.179265e+02 | 0.00000 | 1.800000e+03 | 3.000000e+00 | 0 |
| Alcohol | 4.602861e+00 | 4.052413e+00 | 0.01000 | 1.787000e+01 | 3.755000e+00 | 194 |
| percentage.expenditure | 7.382513e+02 | 1.987915e+03 | 0.00000 | 1.947991e+04 | 6.491291e+01 | 0 |
| Hepatitis.B | 8.094046e+01 | 2.507002e+01 | 1.00000 | 9.900000e+01 | 9.200000e+01 | 553 |
| Measles | 2.419592e+03 | 1.146727e+04 | 0.00000 | 2.121830e+05 | 1.700000e+01 | 0 |
| BMI | 3.832125e+01 | 2.004403e+01 | 1.00000 | 8.730000e+01 | 4.350000e+01 | 34 |
| under.five.deaths | 4.203574e+01 | 1.604455e+02 | 0.00000 | 2.500000e+03 | 4.000000e+00 | 0 |
| Polio | 8.255019e+01 | 2.342805e+01 | 3.00000 | 9.900000e+01 | 9.300000e+01 | 19 |
| Total.expenditure | 5.938190e+00 | 2.498320e+00 | 0.37000 | 1.760000e+01 | 5.755000e+00 | 226 |
| Diphtheria | 8.232408e+01 | 2.371691e+01 | 2.00000 | 9.900000e+01 | 9.300000e+01 | 19 |
| HIV.AIDS | 1.742104e+00 | 5.077784e+00 | 0.10000 | 5.060000e+01 | 1.000000e-01 | 0 |
| GDP | 7.483158e+03 | 1.427017e+04 | 1.68135 | 1.191727e+05 | 1.766948e+03 | 448 |
| Population | 1.275338e+07 | 6.101210e+07 | 34.00000 | 1.293859e+09 | 1.386542e+06 | 652 |
| thinness..1.19.years | 4.839704e+00 | 4.420195e+00 | 0.10000 | 2.770000e+01 | 3.300000e+00 | 34 |
| thinness.5.9.years | 4.870317e+00 | 4.508882e+00 | 0.10000 | 2.860000e+01 | 3.300000e+00 | 34 |
| Income.composition.of.resources | 6.275511e-01 | 2.109036e-01 | 0.00000 | 9.480000e-01 | 6.770000e-01 | 167 |
| Schooling | 1.199279e+01 | 3.358920e+00 | 0.00000 | 2.070000e+01 | 1.230000e+01 | 163 |
Zadbanie o to, aby wykresy były niesamowicie eleganckie:
options(repr.plot.width=24, repr.plot.height=16)
Tworzenie eleganckich wykresów rozkładu zmiennych: warto zwrócić uwagę na wykorzystanie density-plotu dla danych ciągłych i barplota dla danych dyskretnych. Uwzględniono - dla kompletności - rozkład krajów i lat.
for (x in colnames(data)){
if (is.numeric(data[[x]]) && x!='Year')
print(data %>%
filter(!is.na(data[[x]])) %>%
ggplot(aes_string(x=x)) + geom_density(fill="#69b3a2", color="#e9ecef", alpha=0.8))
else
print(data %>%
filter(!is.na(data[[x]])) %>%
ggplot(aes_string(x=x)) + geom_bar() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 6, hjust = 1)))
}
Dodatkowe wykresy prezentujący korelacje pomiędzy Life.expectancy a pozostałymi zmiennymi: pozwalają one zauważyć, że pewne wartości mogły nie zostać poprawnie udokumentowane.
for (x in colnames(data)){
if (x == 'Life.expectancy')
next
plot_data <- data %>% filter(!is.na(data[[x]]) & !is.na(data[['Life.expectancy']]))
plot_itself <- plot_data %>% ggplot(aes_string(x=x, y='Life.expectancy')) + geom_point()
if (x == 'Country')
plot_itself <- plot_itself + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 6, hjust = 1))
print(plot_itself)
}
Z wykresów tych można wywnioskować, że część obsserwacji używała innej skali dla pewnych zmiennych niż wymagana: np. BMI rzędu poniżej 7 (które występowało w krajach o relatywnie wysokiej długości życia) oznacza wagę rzędu ~10 kg dla człowieka o wzroście 1.70 m - można powątpiewać w zasadność tych liczb. Podobnie ma się rzecz z wyszczepialnością; W kolumnie Diphtheria jest bardzo widoczna linia na wysokości ok. x=10; istnieją co najmniej 2 kraje, dla których wartość ‘Diphtheria’ “skakała” pomiędzy wartościami z przedziału 8-9 i 80-90 - Norwegia i Mołdawia.
knit_print(data %>% filter(Country == 'Republic of Moldova' | Country == 'Norway') %>% select(c(Year, Country, Diphtheria, BMI)) %>% arrange(Country, Year))
| Year | Country | Diphtheria | BMI |
|---|---|---|---|
| 2000 | Norway | 9 | 53.3 |
| 2001 | Norway | 91 | 54.0 |
| 2002 | Norway | 93 | 54.6 |
| 2003 | Norway | 92 | 55.3 |
| 2004 | Norway | 92 | 55.9 |
| 2005 | Norway | 91 | 56.5 |
| 2006 | Norway | 94 | 57.0 |
| 2007 | Norway | 93 | 57.5 |
| 2008 | Norway | 94 | 58.0 |
| 2009 | Norway | 94 | 58.5 |
| 2010 | Norway | 93 | 58.9 |
| 2011 | Norway | 94 | 59.4 |
| 2012 | Norway | 95 | 59.8 |
| 2013 | Norway | 94 | 6.3 |
| 2014 | Norway | 93 | 6.8 |
| 2015 | Norway | 95 | 61.2 |
| 2000 | Republic of Moldova | 95 | 46.5 |
| 2001 | Republic of Moldova | 97 | 46.8 |
| 2002 | Republic of Moldova | 97 | 47.2 |
| 2003 | Republic of Moldova | 98 | 47.6 |
| 2004 | Republic of Moldova | 98 | 47.9 |
| 2005 | Republic of Moldova | 98 | 48.3 |
| 2006 | Republic of Moldova | 97 | 48.7 |
| 2007 | Republic of Moldova | 92 | 49.1 |
| 2008 | Republic of Moldova | 9 | 49.5 |
| 2009 | Republic of Moldova | 85 | 49.9 |
| 2010 | Republic of Moldova | 9 | 5.4 |
| 2011 | Republic of Moldova | 93 | 5.9 |
| 2012 | Republic of Moldova | 92 | 51.5 |
| 2013 | Republic of Moldova | 9 | 52.1 |
| 2014 | Republic of Moldova | 9 | 52.7 |
| 2015 | Republic of Moldova | 87 | 53.4 |
Można podejrzewać, że część wyników podano w błędnych skalach - wiedza ta zostanie wykorzystana przy tworzeniu regresora. W przypadku tych dwóch zmiennych można relatywnie łatwo wyodrębnić dane w innej skali; z kolei dla np. Adult.mortality nie dokonano zmian, gdyż linią podziału nie jest stała.
Heatmapa przedstawiająca korelacje między zmiennymi:
cor_matrix <- data %>% select(-c(Country, Status)) %>% cor(use = "complete.obs") %>% round(2)
melt_matrix <- melt(cor_matrix, na.rm = TRUE)
#Wzorzec I: https://r-charts.com/correlation/heat-map-ggplot2/
#Wzorzec II: http://sthda.com/english/wiki/ggplot2-quick-correlation-matrix-heatmap-r-software-and-data-visualization
ggplot(data = melt_matrix, aes(Var2, Var1, fill = value))+
geom_tile()+
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 5, hjust = 1))+
geom_text(aes(label = value), color = 'black', size = 1.5)+
coord_fixed()
Preparacja danych do stworzenia interaktywnego wykresu długości życia w funkcji kraju i roku:
unique_countries <- unique(data$Country)
buttons_countries <- list()
for (i in seq_along(unique_countries)){
buttons_countries[[length(buttons_countries)+1]] <- list(method = "restyle",
args = list("transforms[0].value", unique_countries[i]),
label = unique_countries[i])
}
Wykres długości życia w funkcji kraju i roku:
#Wzorzec: https://stackoverflow.com/questions/63906441/plotly-r-using-colour-and-transformation-with-a-line-plot
p <- data %>% filter(!is.na(Life.expectancy)) %>%
plot_ly(
type = 'scatter',
x = ~Year,
y = ~Life.expectancy,
mode = 'markers',
transforms = list(
list(
type = 'filter',
target = ~Country,
operation = '=',
value = unique_countries[1]
)
)) %>% layout(
updatemenus = list(
list(
type = 'dropdown',
active = 0,
buttons = buttons_countries
)
)
)
p
Preparacja danych dla regresora - w szczególności wykorzystano informację o tym, czy w danej kolumnie zaszedł NaN - tworzono wtedy nową koumnę ‘col_name_is_NaN’ sprawdzającą, czy dana obserwacja ma NaN-a w danej kolumnie, a następnie przypisując medianę istniejących wartości z danej kolumny do NaN-a. Co może się wydawać zaskakujące, to osiągało lepszą skuteczność w regresorze niż wstawienie 0 zamiast NaN-a - może to wynikać z tego, że medianę liczono przed podziałem na zbiór treningowy i testowy. Ponadto wykorzystano obserwację z analizy korelacji między Life.expectancy a pewnymi problematycznie przeskalowanymi zmiennymi.
change <- function(data, name, borderline, multiplier) {
data[[name]][data[[name]] < borderline & !is.na(data[[name]])] <- data[[name]][data[[name]] < borderline & !is.na(data[[name]])]*multiplier
data
}
data <- data %>% filter(!is.na(Life.expectancy))
data <- change(data, 'BMI', 8.5, 10)
data <- change(data, 'Diphtheria', 10, 10)
data <- change(data, 'Polio', 10, 10)
for (x in colnames(data)){
if (data %>% filter(is.na(data[[x]])) %>% nrow == 0)
next
data[[paste(x, '_na', sep='')]] <- rep(0, nrow(data))
data[[paste(x, '_na', sep='')]][is.na(data[[x]])] <- 1
mid = data[[x]]
data[[x]][is.na(data[[x]])] <- median(data[[x]], na.rm=TRUE)
}
Usunięcie powtarzających się kolumn po dodaniu kolumn z NaN-em:
marked = 'Q'
while (marked != '-'){
marked='-'
for (x in colnames(data)){
for (y in colnames(data)){
if (!is.numeric(data[[x]]) || !is.numeric(data[[y]]) || x==y)
next
summa = ifelse (data[[x]] == data[[y]], 1, 0)
if (sum(summa) == nrow(data)){
marked = y
break
}
}
if (marked!='-') break
}
if (marked!='-')
data <- data %>% select(-c(y))
}
## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(y)` instead of `y` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
Wykonanie regresora (warto zauważyć, że regresor nie wykorzystuje kolumny “Country”):
tmp_data= sample.split(data, SplitRatio = 0.3)
train = subset(data,tmp_data==TRUE)
test = subset(data,tmp_data==FALSE)
lm_death <- train %>% select(-c(Country)) %>% lm(formula=Life.expectancy ~ .)
summary(lm_death)
##
## Call:
## lm(formula = Life.expectancy ~ ., data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.8979 -2.1109 0.0173 2.2059 11.8731
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.970e+01 6.681e+01 1.193 0.23324
## Year -1.487e-02 3.333e-02 -0.446 0.65548
## StatusDeveloping -2.127e+00 4.586e-01 -4.638 4.07e-06 ***
## Adult.Mortality -1.780e-02 1.417e-03 -12.564 < 2e-16 ***
## infant.deaths 7.543e-02 1.568e-02 4.811 1.78e-06 ***
## Alcohol 5.049e-02 4.438e-02 1.138 0.25564
## percentage.expenditure 2.256e-04 1.457e-04 1.548 0.12195
## Hepatitis.B -2.800e-02 6.520e-03 -4.294 1.95e-05 ***
## Measles -5.050e-07 1.520e-05 -0.033 0.97351
## BMI 9.851e-02 1.330e-02 7.405 3.16e-13 ***
## under.five.deaths -5.596e-02 1.139e-02 -4.913 1.08e-06 ***
## Polio 1.026e-01 3.251e-02 3.156 0.00166 **
## Total.expenditure -9.412e-03 5.698e-02 -0.165 0.86884
## Diphtheria 4.490e-02 3.236e-02 1.387 0.16568
## HIV.AIDS -4.335e-01 3.166e-02 -13.692 < 2e-16 ***
## GDP 5.362e-06 2.260e-05 0.237 0.81251
## Population -2.231e-09 2.523e-09 -0.884 0.37698
## thinness..1.19.years 2.399e-03 7.893e-02 0.030 0.97576
## thinness.5.9.years 7.808e-02 7.720e-02 1.011 0.31214
## Income.composition.of.resources 5.386e+00 1.021e+00 5.273 1.70e-07 ***
## Schooling 4.919e-01 8.192e-02 6.005 2.83e-09 ***
## Alcohol_na 2.760e+00 1.420e+00 1.944 0.05225 .
## Hepatitis.B_na 7.275e-01 4.204e-01 1.730 0.08391 .
## BMI_na -3.485e+00 1.307e+00 -2.666 0.00783 **
## Polio_na 2.835e+00 2.085e+00 1.360 0.17423
## Total.expenditure_na -1.277e+00 1.295e+00 -0.986 0.32453
## GDP_na -9.976e-01 6.503e-01 -1.534 0.12537
## Population_na 1.261e+00 5.424e-01 2.325 0.02030 *
## Income.composition.of.resources_na -1.116e+00 7.520e-01 -1.483 0.13836
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.661 on 848 degrees of freedom
## Multiple R-squared: 0.8592, Adjusted R-squared: 0.8546
## F-statistic: 184.8 on 28 and 848 DF, p-value: < 2.2e-16
pred <- predict(lm_death, test)
Można zauważyć, że - z pominięciem kraju - najlepszymi predyktorami zmiennej decyzyjnej były zmienne: Schooling, Hepatitis.B, Status, Adult.Mortality, Diphtheria, HIV.AIDS i Income.composition.of.resources - są to też zmienne, które zostały wskazane jako dobrze skorelowane ze zmienną decyzyjną w pokazanej heatmapie. Zapewne nie można argumentować o wpływie tych zmiennych na oczekiwaną długość życia, ponieważ korelacja może informować jedynie o związku, ale nie o przyczynowości.
Wartość \(rmse\):
sqrt(mean((pred - test$Life.expectancy)^2))
## [1] 3.846439
Po przeskalowaniu problematycznych wartości BMI i części szczepień wartości \(rmse\) spadły z 4.11 do 3.8. Przy wykorzystaniu kraju jako predyktora spadłyby jeszcze bardziej (poniżej 2), natomiast jego wykorzystanie można potraktować jako problematyczne; w praktyce można traktować kraj jako pewne ID i podzielić dane kraj po kraju do danych treningowych albo testowych - wtedy predyktor bardzo silnie zależny od kraju kompletnie by sobie nie radził.